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Abstract 

We develop a theory of laser beam generation of shift currents in non-centrosymmetric 
semiconductors. The currents originate when the excited electrons transfer between different 
bands or scatter inside these bands, and asymmetrically shift their centers of mass in elemen- 
tary cells. Quantum kinetic equations for hot-carrier distributions and expressions for the 
induced currents are derived by nonequilibrium Green functions. In applications, we simplify 
the approach to the Boltzmann limit and use it to model laser-excited GaAs in the presence 
of LO phonon scattering. The shift currents are calculated in a steady-state regime. 



1 Introduction 



Photovoltaic phenomena in semiconductors can originate in the built-in or induced 
asymmetry or inhomogeneity of their crystal structures [|l|] . In non-centrosymmetric 
(NCS) crystals, different generation rates for carriers at momenta ±k can be induced 
by asymmetric electron-hole scattering and other processes. The resulting momen- 
tum imbalance generates the so called ballistic current. Recently || , this momen- 
tum asymmetry of carrier generation in semiconductors was achieved by mixing one- 
and two-photon transitions at frequencies 2u Q and uo, respectively. In Ref. [Q, we 
describe the effect in GaAs in the presence of scattering on LO phonons. We have 
also suggested that the current induced by this two-beam coherent control could 
drive intercalated atoms in carbon nanotubes ||. 

In bulk NCS semiconductors, light induced interband transitions of electrons 
in reciprocal space are accompanied by their asymmetrical shifts in the real space 
between atoms in elementary cells. Similar shifts occur if the transitions are induced 
by scattering or in recombination. First realistic description of analogous effects in 
magnetic materials was given by Luttinger ||. The carrier shifts generate the so 
called shift current |7], [9], lQfl , which has in general three components, excitation 
J e , scattering J s and recombination J r , named after the corresponding processes. 
The carrier transitions and the related currents are symbolically sketched in Fig.[l|a. 
The ballistic current and the excitation part of the shift current J e in GaAs were 



investigated analytically in Ref. [11]; J e induced by transitions between light and 



heavy hole bands was also studied ||12|| . A summary of older results about the shift 
current is presented in Ref.[0, but the derivations of the relevant formulas is less 
clearly documented. 

Here, we develop a quantum kinetic theory for the shift current induced by 
laser excitation and scattering. The expressions for the current components and 
the transport equations for carrier populations are derived by nonequilibrium Green 
functions [13| (NGF). In applications, we simplify the approach to the Boltzmann 
limit and use it to study optically excited GaAs in the presence of scattering by 
LO phonons. In our modeling, we consider steady-state excitations by a linearly 
polarized light, and find that J e is reasonably large, while J s is two orders smaller 
and J r is negligible. Therefore, continuous electron pumping through the crystal 
can be achieved. The ultrafast response of J e , without additional saturation and 
relaxation tails from J s and J r , could be useful in opto-electrical applications. 

The paper is organized as follows. In Sec. 2 we describe our model system. 
Section 3 is devoted to derivation of expressions for the current components. In 
Sec. 4 these expressions are further approximated. Numerical results for the hot- 
electron populations and the induced currents in GaAs are presented in Sec. 5. 
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Figure 1: Scheme of the shift current generation in a NCS semiconductor is shown in the inset a). 
The current has three components, related with the process of electron excitation J e , scattering 
J s and recombination J r . They result as a consequence of real space shifts of electrons in the 
elementary cells that undergo the corresponding transitions in reciprocal space. In the inset b) 
we show a small section of a zinc-blend structure in GaAs. The excitation conditions and the 
generated shift currents J e , s ,r are described in the text. 

2 The system studied 

The physics of the phenomenon can be understood from Fig.[I|b, where a segment of 
the zinc-blend structure for GaAs is shown. The valence band states are predomi- 
nantly localized around the As atoms, while the conduction band states are shifted 
toward the Ga atoms. Therefore, if the light is polarized along the (1,1,0) direction, 
the excited electrons transfer from the As atoms at the bottom to the Ga atoms 
in the middle, giving the excitation current J e in the (0,0,-1) direction (negative 
charge). If the light is polarized in the (1,-1,0) direction, electron transitions along 
bonds orthogonal to the light polarization, i. e. from the As atoms at the top to the 
Ga atoms in the middle, generate J e in the (0,0,1) direction. During relaxation, the 
excited electrons (holes) slightly move their centers of charge and stay close to the 
Ga (As) atoms. Therefore, J s is rather small, as we show in our calculations. On 
the other hand, scattering redistributes the carrier momenta, so that electrons at 
Ga atoms recombine symmetrically with holes at all neighbor As atoms, which gives 
negligible J r . 



2.1 The model Hamiltonian 

The space shifts of carriers and the related currents can be calculated from a combi- 
nation of interband and intraband transitions. The length gauge with the elements 
of the position operator x evaluated as in Blount's work |Po| gives a convenient 
basis for the description [|ToT|. Therefore, we model the photo-excited bulk GaAs by 
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the following Hamiltonian (16 



a;k q 

- i e E(t) ■ { ^ E ( G l,k - a«,k) - * E M k ) a l,k «/3,k} 

+ E M «/3( k 7 k - q) a L,k «/3,k-q Oq + fcl 



a;k c,/3;k 

-q/ > 

Q,/3;k,q 



where coupling to LO phonons is added. Here the creation (annihilation) operators 
a a k (a ajk ) describe electrons with the band index a at wave vector k in the Brillouin 
zone. The operator fe q (6 q ) creates (annihilates) phonons with the wave vector q. 
The electric field is E(t) = E +£Jo (t) e _<w °* + E_ wo (t) e +iw °*, where E +£Jo (t) = EI Wo (t) 
are complex envelope functions. Spins are included in the current by a factor of 2 
(see Eqn.(0)). 

2.2 The matrix elements 

The matrix elements x Q!( g(k, k ) for the position operator are defined as JIIJ 



x a/3 (k, k') =iS al3 ^ 5{k - k') + <f(k - k') ^(k) , (2) 
where the functions 

f d 
W( k ) = CpM = J dx < k (x) i — u^x) (3) 

u.c. 

are integrals over the unit cell of the fast components u Q k(x) in the Bloch wave 
functions ^ak(x) = e lk ' x -Ua, k (x). In the following, we use the name r a p(k) = £ a( g(k) 
for the band off-diagonal elements a ^ [3. They are related to the interband velocity 
elements by 

v Q/3 (k) = (a,k\(-ih/m e )V\(3,k) = ^ (k) ~ ea(k) r a/3 (k) . (4) 



The electron-phonon matrix elements [17, 18 



M Q/3 (k, k - q) = M(q) j a p(k, k - q) , 

7 Q/3 (k,k-q) = J dx < k (x) M /3k _ q (x) (5) 

u.c. 

include the structure factors 7 a/ a(k, k — q), which separately depend on the "in- 
coming" and "outgoing" momenta 7 aj a(k, k — q) ^ /(q)- Therefore, the factor 
M(q) in the expression (|5|) gives the rate of electron scattering, while the real space 
shifts of electrons are solely controlled by the lattice M, ITH] , determining the factors 
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7 a/ g(k, k — q) (see Appendix B). In the approximation of a constant LO phonon 
energy huo q m hujQ, the element M 2 (q) is given by [18 



M 2 (q) =M Ml = 2^hu Q (---)■ (6) 

|q| \ £ oo eo/ 

The parameters relevant for GaAs are used; Uuiq = 36 meV, eo = 12.5 and eoo = 10.9. 
3 Description of the system 



We describe the problem by nonequilibrium Green functions jL9[ [20) in the matrix 
form G a p, similarly as in Ref.pj. The shift current can be expressed through the 
interband elements for NGF, which must be known on a semi-analytical level, so 
that the space shifts can be obtained first by algebraic operations. In older works, 
the manipulations have been carried out in the density matrix formalism, but the 
intermediate steps were less clearly presented. More recently, NGF were used to 
investigate the shift current from electron hopping in superlattices [p2TJ , and in other 
photovoltaic phenomena p2l p3[. 



Here, we develop an approach similar to Ref . [pT|[ , and apply it to our system 
in a steady-state regime. The interband functions G a p (a ^ f3) are expressed from 
the differential version of the Kadanoff-Baym equations (19| (KBE) in terms of the 
intraband G aa from the scattering integrals of these KBE. Then G a p are substituted 
in the formula for the shift current, where the mentioned reordering is performed. 
Finally, G aa are calculated numerically from the integral version of KBE M . We have 
not been able to perform this reordering, at least for J s , when G a p was expressed 
from the integral KBE. 

3.1 The differential KBE 



The differential version of KBE for the matrix correlation function G < is [[191 , |2C 
(x = (k, t), integration over \ is implied) 

G i? '°(x,x)V 1 G < (x,x') = (E(x,*) G(ZX)Y , (7) 



where (G^' ) 1 =(G'f^ ) 1 5 a /3 is the inverted free propagator fl9|| . The lesser op- 
eration < can be applied to the functions on the r.h.s of (0) with the help of the 
Langreth-Wilkins (LW) rules fZ§ (AB)< = A<B A + A R B<. The KBE can be also 
written in a form, where (G A '°)~ acts on G < from the right side, and, in the 
scattering term, X is interchanged with G. 

Subtraction of these two KBE gives @ (T = (t + t')/2) 

ih dG<a ^ X) - (e Q (k) - e ,(k')) G^(X,X) 
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= S /ia? (x, X) G^{x, x) ~ G^(x, X) Efflux, x) 

+ (e«;«7(x,x) G^fax')) - ( G a7(x,x) S s; ^(x,x')) • (8) 

In the scattering integrals on the r.h.s., field £/ ;a /3 an d scattering E s;a/ 3 self-energy 
functions are introduced Q. The terms do not commute because of the time, mo- 
mentum and band indices. The momentum and band non-commutativity ultimately 
leads to the shift currents. The time issue should be less serious, if pulse fields with 
slow envelope functions are applied to the system, where gradient expansions in 
terms of the difference time r could be performed around the CMS time coordinate 



T [ 19|, psfl . We adopt this approach here and include only the zero order terms. 



3.2 The self-energy functions 

In order to obtain the field self-energy £/ ;a /3 for the Hamiltonian (|l|) in the length 
gauge, it is necessary to perturbatively expand the Green functions in terms of the 
light excitation, similarly as in Ref.Q]. Direct calculation gives 

E /;Q/3 (k,k',t,t') = -% e 8(t - t) 5(k - k') 

which has zero correlation parts Ey.'^(k, k , t, t ) = 0. This self-energy acts as an 
operator, due to the presence of derivatives. When surrounded by functions with 
two momenta variables, it differentiates the front (back) function over k (k'), and 
set k = k . The derivatives can be shifted from one side to the other side by a 
partial integration, which changes their signs and combines the two with a prefactor 
1. We can also introduce its steady-state form 

£± Q/3 (k,k') = -ieS(k-k')E ±a;o .{i (^"D^-'Wk)}, (10) 

which can be obtained by a Fourier transform of the Dyson equation in the frequency 
domain ||19|| . The frequency argument of the Green functions following Sy. a/3 are 



shifted by =Ftu (see Ref.0). 

Finally, we have to write down the correlation function for the electron-phonon 
self-energy, which we use in the self-consistent Born approximation 



£< Q/3 (k, k', t, t ) = M Q7 (k, k - q) M 5/3 (k' - q, k 



x G< 5 (k-q,k'-q,M') D<(q,t,t) . (11) 



Here, the non-locality in momentum and time as well as all kinds of interband 
transitions are included (summation over q and repeated band indices). 



G 



3.3 The total shift current density J 



The shift current density J can be expressed in terms of the interband velocity 
elements Vp a (k) and the k-diagonal elements of the correlation functions G<p(k, k). 
If the perturbing electric field is homogeneous in real space, i.e. invariant with 
respect to the lattice translation of the crystal, then we would naturally expect 
that these functions are k-diagonal G^ik, k') = G^ /3 (k) 5(k — k') (2ir) 3 . In reality, 
Eqn.(||) for the Hamiltonian (jl]) in the length gauge should be further transformed 
|20H , to explicitly give zero off-diagonal elements, and this could generate additional 



terms. In our problem these transforms would be complex, so that we make the 
Ansatz that Eqn.(||) gives momentum diagonal G^p- 

The steady-state shift current density J can be expressed as follows 

dhu r dk 



26 J J J^Y S ^(k) G<p(k, U ) = J e + Js + Jr , (12) 



where the prefactor 2 encounters the spins. The components J e ,s,r, which represent 
contributions of the three processes depicted in Fig.|l|a, can be obtained, if vp a from 
(|J) and the solution G<p of Eqn.(H) are substituted in the expression (|I2"D. In the 



steady state, the first term on the l.h.s. of Eqn.@ is absent. The remaining G^p, 
expressed through the scaterring integrals on the r.h.s., can be substituted in (|12"D, 



which was also used in Ref . pi] . In transient situations, the time derivative of G^p in 



Eqn. (R) must be also included, but we are postponing this problem to future studies. 
Terms with equal band indices a = [3 in (|T2| ) would give the ballistic current density, 
calculated in Ref.|Q for excitation by two laser beams. 

3.4 The excitation current density J e 

Let us first find the expression for J e from the two terms with on the r.h.s. 

of Eqn. (H) . Since J e results in the second order of the laser field, ^f- a p in those 
terms must be combined with the first order interband correlation functions Gf. Q/3 ~ 
(G . aa G .ppj , expressed through the interacting Green functions in the 

absence of laser field 0] G 0]aa . The resulting J e has the form 

r duj r dk ^ 

x ( s / ; a 7 ( k ; ki) (G ;7 7 (ki, k 2 , u) Sj ;7/3 (k 2 , k 3 ) G .pp(k 3 , k,u± u ] 
- (G 0; aa(k,ki,a;) Ej. a7 (ki,k 2 ) G 0;77 (k 2 , k 3 , u ± ^ )) < £/. 7/ g(k 3 , k) \ ,(13) 

where we neglect vertex corrections to Gf. af3 , which can contribute especially in 
transient situations ||, and assume that G . a p ~ 0. We also write two momenta 



< 
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in the equilibrium Green functions Go. aa (k, k'), in order to perform the deriva- 
tives from the field self-energies £j ;Q/3 (k, k'), even though they fulfill G - aa (k, k') = 
G . aa (k) 5(k — k') (27r) 3 without the above Ansatz. In Appendix A, we show [TJ| 
that after algebraic manipulations only elements with 7 = (7 = a) remain in the 
first (second) expression of Eqn.flTBD. 

In total, we obtain the c-component of the vector for the current density J e 

G ; aa (k,Uj) G 0]( lf3(k,UJ ± UJ )j , (14) 

where sum over a and b components is performed. The sign convention and the 
terms ^ a . c (k) are described in Appendix A. For weak scattering, the expression 
([14]) agrees with density matrix calculations JIEfl . In this formalism also a nonzero 



virtual term below the band gap can be traced in an expression analogous to (|H|). 
It can be canceled by a similar term, with opposite sign, resulting from the band 
diagonal contribution v aa G< a to the current. 

For a light linearly polarized in the 6-direction, the product of the components 
E±oj w hh r b p a . c (k) r^k) can be further simplified. If the matrix elements are used 
in the form || [TO] r b ^ a (k) = \r b ^ a (k)\ e l ^f a( ' k \ we can arrange J e as follows 

Je - 2 e J 2^ J J^f lE ^\ R e ]P M Mk)| + - dk ) 

x (G 0;aa (k,u) G .pp(\t,uj ±u )) , (15) 



where we have introduced the excitation shift vector R*.^ with the c-component 

^a(k) = %^ + C(k) - ^(k) • (16) 

From the expression (|16|), it follows that the vector is invariant under the phase trans- 
formation H [10], |TJJ -?/>/3k( x ) — > e* 9n ^ k V/3k( x ) of the Bloch functions ^gk(x), where 
6 n (k) are arbitrary well-behaved functions. It is also ant i- symmetric Rg.^ Q (k) = 
— Rg. a/3 (k) in the band index a «-> /3, since r a/3 (k) = (r^k))*. Note in (|T5|), that 
the anti-symmetric Rg. a/3 (k) combines with the imaginary (spectral) values from 
the Green function product, while the band symmetric derivative 9|r^g(k)| 2 /<9k 
combines with its real (main) part, representing renormalization effects due to in- 
teractions. An analogous situation occurs in the scattering shift vector R s;/ g j a(k, k') 
(see Appendix B) and in the scattering current density J s , discussed below. 

3.5 The scattering current density J s 

The formula for J s can be obtained, if the rest two terms from the r.h.s. of Eqn.(|[), 
with the electron-phonon self-energy in ([IT]) , are substituted in the expression 
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In this work, we consider only intraband relaxation, so that band-diagonal Green 
functions are used in J s . We also assume that they depend on one k variable, in 
accordance with our Ansatz. 

Let us concentrate first on the two terms r (E R G < — G i? S < ), in the expression 
resulting from (12). They can be written as follows 

£ r/3Q (k) (x« p (k,u) G< p (k,u)-G* a (k,u) X< p (k,u)) 

a^P V 7 

= i>(k)M Q7 (k, k - q)G* (k - q, u - u)D>(q, to - cu)M 7/3 (k - q, k)G^(k, u) 

- 7 r/3a (k)M a7 (k, k - q)G< 7 (k - q, u - u)D R (q, u - Qo)M lP {k - q, k)G^(k, u) 

a^P; 7 

- £r^(k)G*(k,u;)M^(k,k-q)^ , 

(17) 

where the band arguments have been shifted in the last term. In the propagator 
part of the electron-phonon self-energy, the formula H R oc (GD) R = G R D > — G K D R 



has been used, which is valid for functions with this argument ordering [24, |25 



The first term on the r.h.s. of (17) , resulting from G R Z) < , looks similar to the last 



term, but the two differ in the type of phonon correlation functions D> , D K and 
the energy-momentum arguments in Gpp, G R y and M 1 p. 

This can changed, if the substitution k = k — q and u = u — (o, followed by 
uj" = —id, is used in the last term of (ff7|) ; we take advantage of the fact that the Bose- 
Einstein distribution satisfies n BE {uj) = 1 / (exp(huj / kT) — 1) = —(1 + n BE (— uj)). 
The sum of the first and last terms is 

i>(k) M Q7 (k, k - q) - Y M /3«( k ' k - q) r «7( k - q)) M 7/s( k - q. k ) 

a^P; 7 7^Q; P 

x (k - q, u - Co) D>(q, u - u) G< p (k, u) . (18) 

Here the band indices at the electron propagator and correlation function can be set 
equal 7 = /3, because intraband relaxation is considered. The commutator in the 
large bracket, multiplied by the element M^pik — q, k), can be rewritten in terms of 
the scattering shift vector R^^k, k — q) and a renormalization correction, as shown 
in Eqn.( p73|) . If interband relaxation is also considered, the resulting formulas can 
be generalized analogously. 

The other two terms r (£<G A - G<S A ) in flTJ) can be reordered similarly, if 



only the term G A D < in S A oc (GD) A = G A D > - G K D A is again considered. The 
remaining terms — G < D R , — G < D A from {GD) R , (GD) A , respectively, can also be 
summed to this form. Here, the imaginary part of the phonon propagator appears, 
instead of the electron propagator. Its real part would give further renormalization 
contributions, if phonon scattering was considered. 
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In total, the scattering shift current density J s in steady state is equal to 

r dw f dk ^ 
= " 4 e / K= LiE 



. 2tt J (2vr)3 ( 

x jk if (k,k-q) |M^(k,k-q)| 2 ImG^k-q.u-w) 

l/ aiM^k-q.k)! 2 <9|M A8 (k-q,k)| 2 X \ 

" 2I dk + — dik^) — J Re ^( k -q^-^)J 

x £>>(q>-u;) G|^(k,c;) + R s; ^(k,k-q) |M^(k,k-q)| 2 

x l(G^(k-q,u-u)-G^(k-^,u + u))G^(k,u)\. (19) 

The intraband correlation functions Ggg can be obtained from the solution of the 
KBE similarly as in Ref.Q. If the hot-electron population is small, then it is possible 
to keep in fli~9|) only the terms linear in Ggg. 

The renormalization terms in ([19]), with derivatives of matrix elements, give 
non-classical corrections to the shift current, due to quasiparticle broadening of the 
carrier spectra. Note however, that J s is practically independent on the scattering 
rate, at least for weak scattering, since R s;/ 3/3 is determined by the crystal structure 
(see Appendix B), and the relaxation carrier flow is equal to the injected flow. 
From the same reason, the temperature dependence of J s is also weak. Expression 
analogous to Eqn.([L4|), ( |19"D can be derived for the recombination current density 
J r , but in Sec.V we argue that J r is usually negligible. 

4 Approximate solution for the shift currents 

Here we approximate the expressions for J e and J s to the Boltzmann limit, which can 
give reasonable results in materials with a moderate scattering. First we simplify J e 
in flT3|), where the light orientation is chosen along the (1,1,0) direction, appropriate 
in GaAs. Next, we approximate the equations for the correlation functions Ggp, and 
insert their solution in Eqn. (|19|) for J s , together with the shift vector R s -pp in (B.3|) . 

4.1 The excitation current density J e 



In numerical studies, it is convenient to evaluate J e from the expression (|14"D even 
for a linearly polarized light, so that ab-initio calculations of the elements rg Q (k) 
can be performed along the usual crystal axes. Once evaluated, r / g Q (k) can be 
used to obtain the elements r^^k) from the formulas (3.36-3.37) in Ref . |15| . At 
weak scattering, renormalization effects contained in the derivatives d\r^(k)\ 2 /dk 



from Eqn.flI5p can be neglected, which, in Eqn. (13), is equivalent to taking only 
the imaginary part of r^ Q . c (k) r^Jk), anti-symmetric in the band indices. Then, 

G< v G^ c (a = v, P = c) there can be combined with — Gf c G< v (a — c, (3 — v), 
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which gives G< v (—2 i Im G^ c ) = i A vv A cc , where A aa are the electron spectral 
functions (§]. The terms with exchanged band indices are zero, since in equilib- 
rium the conduction band is empty; i.e. Gf c (k,u) = A cc (k,u) npniyj) « 0, where 
ufd(^) = l/(exp(hu/kT) + 1) is the Fermi-Dirac distribution. 

The k and cj-integration in the two spectral functions A vv A cc is performed as in 
Ref. Q. A parabolic approximation for the bands is considered in these integrations 
(not in r ay g(k)), with wave vectors tuned to the resonant values k™ es , used to describe 
electrons that relax between the energy levels E™ es = E® es + uTuoq (laser tuning to 
E® es ). The notation for angles is related to the electric field E^ = E wo / \/2 (1, 1, 0); 
the angle G (0, 7r) is between the direction (1,1,0) and the k™ es -vectors, and the 
angle 9 G (0, 2ir) is in the plane orthogonal to the (1,1,0) direction. The nonzero 
^-component of the excitation current density is (77 = (k® es , 0, 9)) 

jz _ ■ ^res c I rp |2 

x ^ [ g sin(0) (r^) rlM + r^M r* vc (rj)) ■ (20) 

where \l cv = m c m v /(m c + m v ) is the effective electron-hole mass, and the factor in 
the bracket is purely imaginary. If the light polarization is in the (1,-1,0) direction, 
the signs at both terms with x arguments change. Since the integral is the same, 
the current is opposite, in agreement with Fig.[l|b. 

4.2 The scattering current density J s 

In order to approximate J s in Eqn. ([l9D , we need first the transport equations for 
the intraband correlation functions GqJIl, u). Here we describe the hot-electron 
relaxation by the integral KBE derived in Ref.0]. In the weak scattering limit, they 
can be approximated to the integral Boltzmann equation (IBE) ||]. For simplicity, 
we also describe the hot carrier population only in the conduction band. 

A laser light polarized along the (1,0,0) direction, used in Ref.[§|, gives a nearly 
isotropic hot-electron distribution in the plane orthogonal to this direction, which 
is sufficiently described by the angle <fi. In the present work, polarization along 
the (1, 1, 0) direction gives a population squeezed in the (0, 0, 1) direction by about 
25 % (see FigFJ). Moreover, the shift vectors are also very anisotropic (see Fig.|3|b), 
so that both angular variables 0, 9 must be used. The field self-energy for interband 
transitions along the (1,1,0) direction can be obtained from the expression (|10"D in 
the form 



S+ CT (0, 9) = -e (r x cv (k° res , 0, 9) + r» 0, 0)) Ej^2 . (21) 

In the scattering self-energy S s;cc , the matrix elements M cc (k, k') use the linearized 
structure factors 7 cc (k, k') from (|B.10j ). Since intraband relaxation is considered, the 
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square of the exponential prefactor from ( |B.10| ) gives 7 cc (k, k') 1, i.e. M cc (k, k ^ 
M(k-k') from (§. 

Then using the steps in Ref. M, we arrive at the steady-state IBE 

o p r (n) ( 

fccM,e) = y ' \ |s+ OT (0^)| 2 ^ 6, 



cv u n0 



™c „ f 2w d6 f n dd . . -. ( 1 „ , , - -, . . 

+ T M « I Sis s '"« ( iift.-ts-p W"" 1 .*.') »'M 

+ ,,„ ^Tjj /«(« + 1,0,9) (1 + %W)) , (22) 

res res I ' ) 

where the distribution function on the n-th level (e* = E™ es ± nhojQ/2) 

Un, 0, 6) = f ^ /°° ^ fc 2 G< cc (fc, 0, (9, W ) (23) 

Jen Z7T JO 27T 

is defined from the second order G K ~ G^/2!, expanded in terms of the interband 
field self-energy S/ ;cy ; the prefactor 2 in Eqn.(p2|) cancels this 2!. Here r D (n) = 
—%/(2 ImE^, s (n)) is the particle relaxation time. It is still necessary to add in 
Eqn. (|22|) radiative transfers of carriers between the bands, but this is practically 
not reflected in J e , and J r is also negligible. 

To calculate J s , we also need a simpler expression for the scattering shift 
vector R^^k, k') in ( |B.4| ) . A tractable approximation ||, |10[ can be obtained if 
this is linearized in terms of the wave vector difference k — k , around the center 
ko = (k' + k)/2. For large differences, where the errors can grow, the contributions 
to scattering are small, since the elements |M(k — k')| 2 decay as |k — k'|~ 2 . Direct 
algebraic manipulation gives the scattering shift vector in ( |J3. 1 1| - |BTT2D . In the current 
formula (|T9|), we expand the vector H s -pp in each considered wave vector k in the 
Brillouin zone. 

Finally, the expression (JE]) for J s can be rewritten in the Boltzmann limit, 
where App(k.,u) = —2 Im Gpp(\t, uS) ~ 2-n5{Tiuj — e / g(k)). Since the correlation 
function -D > (q, uj) for free phonons is also sharp, it can be easily convoluted over 

frequency and momentum with App(k, u). If we consider in (|T9| ) only the term linear 

< 

A3' 



in Gg 3 , and neglect renormalization factors, then the approximate J s cc results 



ft „ JO Z7T JO Z7T 

h.n+1 _ \ 

+ W~ n ; £5 r n+1|2 (Kes - Kts 1 ) x n B (uj Q ) . (24) 

I res res I / 

Here the distribution from (|23|) has been used, where the integration over uj and k = 
|k| from (19) is already performed. A more consistent approximation in Eqn.(^4j) can 
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be obtained by evaluating the vectors Q c at the points (k™ es + k^)/2. Numerically, 
it is easier to average fi c , evaluated at the side points k" es , k^. 

5 Numerical results and discussions 

In applications, we consider a typical experimental configuration for GaAs; the laser 
intensity is I = 10 6 W/cm 2 at the light energy Hujq = 2.1 eV (band gap E g = 1.5 
eV), with light polarized along the (1,1,0) direction. We calculate the steady-state J e 
from Eqn. (|20|) , and J s is obtained by solving Eqn. (|22|) for the intraband distribution 
ff3i3(n,<j),9), which is then used in Eqn. ([24]). The recombination current density J r 
is also discussed. 

We consider a model of GaAs with 10 bands, without spin-orbit coupling. The 
matrix elements are calculated ab initio within the density functional theory in 



local density approximation using a plane wave-pseudopotential approach p6[ . They 
are obtained at the resonant momentum k® es (resonant value for the light energy 
Uujq) and at about 40 energy (momentum) levels in the band c, corresponding to 
resonant values of LO phonon processes (E™ es at momenta k™ es ). The elements are 
found on a mesh (0j, 9j), which can be reduced due to symmetry to a size 10 x 20 
points in the region (p = (0, vr/2), 9 = (0, 7r), giving in total around 100 Mb of input 
data. The calculated current density J e corresponds to the excitation between the 
heavy hole bands (v = 3, 4 in Eqn.(^l|)) and the lowest conduction band (c = 5), 
while J s is illustrated only in the band c = 5. 

5.1 The electron distributions 

In Fig.|2], we show cross sections, orthogonal to the (1,1,0) direction, through the 
electron population f cc {n, <fr, 9) in the conduction band, multiplied by sin(0) as in 
Eqn.(|24|). The solid and four dashed lines on each plot correspond to the angles 
<f) = i it/ 10 (i — 5 — 1), while the shapes of the ovals give the ^-distribution. The 
presented levels are n = 0, —1, —2, and the temperature is T = 300 K. Note, that the 
population {n = 0) is squeezed in the (0,0,1) and (1,1,0) directions, where the last 
is seen on the fast decays with of the ovals for n = 0. The population squeezed in 
the plane orthogonal to the (1,1,0) direction, around <fi m 0, contributes more to the 
shift current than the population in the plane orthogonal to the (1,-1,0) direction. 
Therefore, these contributions with opposite signs do not cancel, and .11 can be 
nonzero. For a light polarized in the (1,-1,0) direction, the population is squeezed 
in the orthogonal direction, according to the change in \r x cv ±r^| 2 from (^), and J| 
changes sign. At lower levels both squeezing become relaxed, which visually shrinks 
the size of the population. Relaxation of anisotropy in hot carrier distributions was 
also studied in Ref. Wi 
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Figure 2: Crossections through the population f cc (n,<f>,9), multiplied by sin(</>), in the planes 
orthogonal to the (1,1,0) direction. The values for levels n = 0,-1,-2 show relaxation of the 
squeezing along the (0,0,1) and (1,1,0) directions. The solid and four dashed lines correspond to 
the angles cf) — i 7r/10 (i — 5 — 1). 



5.2 The shift current density J e 

In Fig.|3|a, we show the angular dependence of A rr = Im [r x cv . z r v vc +r% v . z rj c ), mul- 
tiplied by sin(0), from the expression ( |20"D for J e . It is calculated for the angles as 
in Fig.[| and for the light polarized in the (1,1,0) direction. If the light is polarized 
in the (1,-1,0) direction, A rr looks the same, but the prefactor in Eqn.(^) changes 
sign. This reflects the fact that the structure of GaAs allows pumping from different 
atoms in both situations, but it gives J e = for excitation by unpolarized light. For 
the present excitation, we obtain from Eqn. (|20"D the value J* ~ —45 A/cm 2 . This 
agrees in sign with Fig.|I]b (negative charge of electrons) and in value with Ref.JTB 



5.3 The scattering current density J s 

Evaluation of J s is sensitive to numerical errors, resulting from approximate integra- 
tions on the finite mesh for 0, 9. We have corrected three problems in the calculation 
of J s -cc- First, it is the flow between level n and nil, second the conservation of 
homogeneous distributions (used as an initial test) in scattering between level n and 
nil, and third the fact that nonzero shift currents might result even for homo- 
geneous distributions, which is a non-physical consequence of the rough mesh and 
other approximations. 

In Fig.|3|b, we show the ^-component of the linearized scattering shift vector 
R z s . cc from the level n = 0; the cross section is orthogonal to the (1,1,0) direction at 
(ft = 7r/2. The behavior of R z s . cc can be appreciated, if we substitute the difference 
of wave vectors in ( [B.ll[ - |rTT2] ) by — k. This because, in average, the wave vector k 
of the excited electron becomes scattered in the — k direction, even though the size 
k® es of — k is about an order of magnitude larger than the difference k® es — kf^ s . The 
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vector R z . cc changes a sign, if the cross section is orthogonal to the (1,-1,0) direction, 
proving zero current for a homogeneous distribution. Because of the symmetry 
specified above, we can limit calculations to one quarter of the Brillouin zone. Note 
that R z . cc is larger in the horizontal direction, where the squeezed population in 
Fig.|] is also larger, so that J z s . cc becomes increased. 



Arr((p,G) sin(cp) [nm] R?, C c(cp=7i/2,e) [nm] 

-0.02 - 0.23 




Figure 3: Inset a) shows the angular dependence of the factor A rr = Im {r£ v . z r\ c + rf! v . z r% c ), 
multiplied by sin(</>), calculated from the level n = 0. The crossections are in the plane orthogonal 
to the polarization direction of light (1,1,0). It is opposite if light is polarized in the direction (1,- 
1,0). This proves the mechanism in Fig.^jb, with zero current J* for excitation by unpolarized light. 
In the inset b) we present the approximate scattering shift vector i?f. cc (</> = tt/2, 8) w (— k x fl c ) z , 
calculated from the level n = 0. The crossection is as in inset a). For the orthogonal crossection, 
the contributions only change sign. Thus J% cc = for homogeneous electron distribution, resulting 
approximately for excitation by unpolarized light. 

In Fig.[|a, we present contributions to J z s . cc from different electron levels in the 
conduction band, as induced by phonon scattering |4j] . The temperature is T = 300 
K, and the dashed, solid and dash-dotted lines correspond to the light excitation at 
the energies corresponding to the levels n = —2, 0, 2 (used the same injection rate). 
The three results are very different, since the shift vector R z s . cc and, consequently, 
J z . cc changes sign around level n = (see also Fig.|). This effect is related to the 
form of matrix elements, and it appears by chance close to the level n = 0. 

The contributions to the current density J z . cc are summed and presented in 
Fig.^b in the range T = 50 — 300 K. The excitation energies are labeled by the level 
number n = 0,5, 10. In this interval the response changes sign, as expected from 
Fig||a. J z . cc is about 3 orders of magnitude smaller than J z , but it increases for 
larger excitation energies (see Fig||). For excitation around n = 0, J z . cc is nearly 
temperature independent, while away from n = it grows in size with T. In general 
the increase is small, since J z . cc depends weakly on the scattering strength (see 
discussion at Eqn.(|P^)). 

In Fig.||, we show the dependence of the ratio J Z . C J J z on the excitation energy, 
labeled as in Fig.|]b. As already mentioned, the ratio is negative close to the T point, 
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Figure 4: Inset a) shows contributions to the scattering shift current density J*. cc from different 
"electron levels" , due to scattering on LO phonons at T = 300 K. The dashed, solid and dot-dashed 
lines correspond to tuning the laser to the energy of the level n = —2,0,2. The change of sign of 
the contributions is due to the fact that R z s . cc at angles 4> — ir/2 changes sign at these excitation 
energies. In inset b) we present the temperature dependence of the scattering shift current density 
Jg. cc calculated as a sum of its components in inset a). The situations correspond to tuning the 
laser energy on the levels n = 0, 5, 10. J% cc changes sign and slightly increases with temperature, 
especially at larger n. 




n 

Figure 5: Ratio of the scattering and excitation shift current densities J s -cc/Je, calculated as a 
function of excitation energy, labeled in phonon levels. The solid (dashed) lines correspond to 
T = 300 K (T = 50 K). The ratio changes sign around n = from positive to negative and 
increases in value up to 1% for the present excitation energies. 
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it changes sign around n = and approaches 1 % at n = 15. Here we stop our 
calculations, since the departure from the parabolic band approximation is large. 
The solid (dashed) curves correspond to T = 300 K (T = 50 K). Since J* roughly 
increases by 25% in this energy region, J z s . cc is responsible for the increased value 
in the ratio. Note also that at higher excitation energies the current density J z s . cc 
is opposite to J|, which could be intuitively expected. The ratio can be larger for 
relaxation to/from local side minima in the band structure. 

5.4 The recombination current density J r 

Finally, let us shortly discuss the recombination current density J r . Since the mo- 
mentum relaxation time r p is several hundreds fs, and the recombination time r r 
is about 1 ns, only a tinny fraction of carriers recombine before randomizing their 
momenta. The hot carriers first isotropically fill the Brillouin zone, and in particular 
the distribution is practically the same in the (1,1,0) and (1,-1,0) directions. It is 
good to recognize that J e and J r resulting from transitions between two particular 
states have opposite signs. Therefore, contributions to J r from the population in 
the above directions are the same in magnitude, but they point in the (0,0,-1) and 
(0,0,1) directions. This means that J r is practically zero, as suggested in Fig.|I]b, 
where the recombining electrons going from Ga atoms do not distinguish between 
the As atoms. The same observation is expressed in Ref. |IJ in a more general way, 
namely, that J r is zero in non-pyroelectric materials. 

Recombination through trap sites from impurities would give nonzero micro- 
scopic shift currents, but its macroscopic value J r averaged over all impurity posi- 
tions should be negligible. It would be also interesting to study the shift current 
related with transitions to excitonic levels. Since these levels are energetically below 
the band edge, it can be expected that they give different (smaller) shifts in real 
space than in the interband transitions. These shifts would also vary level from 
level, and J r could be also nonzero. 

6 Conclusion 

We have theoretically investigated laser beam generation of the shift current den- 
sity in bulk NCS semiconductors. The excitation, scattering and recombination 
components J e , s ,r reflect asymmetric carrier flows in elementary cells, induced by 
the relevant processes. The system was described by the NGF methods, which, in 
combination with the length gauge, gives a consistent approach to this problem. 
Expressions for J e s were derived for steady-state excitations, in terms of the carrier 
transition rates and the space shift vectors R e , s - 

For practical purposes, we have simplified the formalism to the Boltzmann 
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limit and demonstrated a tractable numerical scheme. Within this approximation, 
we have described optically excited GaAs in the presence of LO phonon scatter- 
ing. Light induced electron transitions between the lowest conduction band and the 
nearest heavy hole bands were considered. For light polarized in the (1,1,0) direc- 
tion, the excitation current density J e in the (0,0,-1) direction was obtained. The 
scattering current density J s;cc was calculated from the hot-electron distribution in 
the conduction band. It is about two orders of magnitude smaller than J e , and the 
recombination current density J r can be neglected, since the momentum relaxation 
causes that electrons recombine with the same strength to atoms placed in opposite 
directions. Therefore, steady-state pumping of electrons across the crystal can be 
realized. 

The shift current has potential applications in ultrafast photo-electric devices, 
since the response is practically instantaneous. Modifications of the shift current 
might be observed in low dimensional structures, heteropolar nanotubes |28| or 
molecular systems. 
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Appendix A 

Let us first compare current expressions from Ref.fl6| with our results. The sum 

T (2) 

of the term (Ji ntra ) and d (P inter ) /dt, there, corresponds to the excitation shift 
current density J e from our Eqn. ([T2|) , if scattering is neglected. The term {Ji ntra } n 
there corresponds to the excitation part of the ballistic current; Eqn. (|T2|) with equal 
bands a — f3, and no transport vertex corrections to G< a . 

Next, we show how Eqn. ([Hp can be obtained from Eqn.flTBD, in analogy to 



Ref . fllql . Two situations can be considered in (|13|) , where either one of the two 
involved field self-energies ^f- a p is interband and the other is intraband, or both 
of them are interband. In the first case, the interband self-energy in Eqn. (|T3D is 
between the two Green functions, since these are from different bands (empty/full), 
to give real transitions. The intraband self-energy in front or back can be switched 
by partial integration (see comment at (|J)) to give a derivative, with a prefactor 1, 
over the front or back variable in the first or second term in Eqn.flTBD with Gf. aj3 . 
These two can be combined into a derivative dGf. a g(]s., k)/9k, and transferred by 
partial integration to the derivative in the prefactor [|TjJ — <9r aj g(k) / dk. The resulting 
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contribution to J e in (p~3|) is 



X ( G , ;aa(k,U;) G ;/3/3(k, W ± CJ ) ) , (A.l) 



where we employ the definition ||T6[ 

r Pa M = - i (Q„00 - C(k)) r^(k) . (A.2) 

In ( [A.l]) the sum over bands a ^ (3 selects real transitions between valence and 
conduction bands, tuned by the laser frequency u . Then, the top (bottom) signs 
correspond to a = c, (3 — v (a = v, (3 — c). 

The second situation adds the expression rg a (k) J2j^a,i3 [£a 7 (k) £ 7j a(k) — 
£L,(k) £ 7/3 (k)] to the bracket (— i r^ a;a (k) ^^(k)) in ( |A.1|) . These terms can be 
largely simplified with the following identity [|l^] 

r^ ;6 (k) - r^. a (k) = i £ [^(k) e»,(k) - C 7 (k) &(k)] , (A.3) 
which can be used to exchange the vector (b) and derivative (a) components in 



r a/3;a(k)- If we substitute ( |A.3|) in (|A.1|) , and shift arguments, then the square 



brackets cancel and the expression fllij) results. 
Appendix B 

Here we derive the expression for the intraband scattering shift vector R^^k, k'), 
used in Eqn.(|T9D, and linearize it in the wave vector difference k — k . 

We assume |3|, |IT]] that the commutator of the operator x in (Q) and M in 
(^) vanishes, i.e, [x, M] = 0. In the Bloch representation, we can easily sum over 
intermediate states and momenta. Upon separating out the terms of the position 
operator that are band-diagonal, we arrive at the following sum rule: 

£ i>(k)M Q7 (k,k')-E M^(k,k')r Q7 (k') 

= ~ l im + w) M/37(k ' k ' } + M/37(k ' k ' } ^ 77(k ' } " ^ (k) ) • (R1) 

Aided with Eqn.( |B.lD and writing the matrix elements Mpp(\t, k') in the form 

%(k,k') = |M^(k,k')| e l ^( k ' k ') , (B.2) 
it is now straightforward to simplify the bracketed terms in Eqn. (|18|) for 7 = (3: 
r^(k) M a/3 (k,k') - £ M Pa (k,k') r^(k')) M^(k',k) 
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= ir^^k') \M^,k)\ 2 'l {dkc + d^) l M ^( k '' k )| 2 • (B-3) 
Here the c-component of the scattering shift vector R^^k, k') is 

KUk, k') = (0 + 0) ^( k , k') + e^(k') - «k) . (b.4) 



It has properties similar to the excitation shift vector R e;a/ g(k) in (16). R s; /3/3(k, k' 



is invariant under phase transforms of Bloch functions, and anti-symmetric upon 
exchanging k and k', i.e. R s; ^(k, k') = — R s;/ g / a(k', k), which follows directly from 
Eqn. (|B.2D . Therefore, the real (first) part in Eqn. ( p73|) is also anti-symmetric, while 
the imaginary (second) part is symmetric, and these two give different contributions 
to the scattering shift current density J s . 

Numerically, it is more convenient to evaluate the shift vector directly from 
from the matrix elements M a/ g(k, k'), which are accessible to ab-initio calculations. 
Then, considerable care has to be exerted in order to maintain the invariance under 
phase transformations. We make use of the specific form of Mpg(k, k'), with the 
structure factors 7/3/3 (k, k') given in Eqn.(|), to evaluate i?^ ;/3/3 (k, k') from 



1 T^(k,k') ( A + JL) 7/3/3 ( k;k ' 



.| 7 w(k,k')i 2 PP ' V^ c M e 
+ ^(k')-^(k). (B.5) 

Since 7/3/3 (k, k') only depends on the band structure, it is clear that the shift vector 
i?s;/3/3 is an intrinsic property of the material and does not depend on the nature of 
the scattering mechanism. 

Eqn.( |B.5l) can be used once 7/3/?(k, k') are found. For practical evaluations it is 
much easier to obtain their linearized form in the difference k — k . We can evaluate 
ii/3k(x) and Ug k '(x) as well as their derivatives from ii/3k (x) at k = (k + k')/2 using 
kp-perturbation theory ||29|| . This allows us to obtain corresponding perturbation 
theoretical expressions for 7 1 g / g(k, k), £0/3 (k) and the shift vector i?^ ;/3/3 (k, k'). 

The linearized «/3k(x) results in the form 

M/3 k0+Ak (x) = e -^( k0 > Ak I u/3 k0 (x) - i Ak • £ r a/ 3(k ) n ak0 (x) I , (B.6) 

V J 

where Ak = k — ko. The exponential prefactor reflects the change of the phase in the 
Bloch wave W/3k +Ak(x), which is otherwise completely undetermined. As mentioned 
above, this phase is of paramount importance, since we have to take derivatives with 
respect to k from the Bloch wave M/3k +Ak( x ) (see Eqn.( |B.5D ). Therefore, we have to 
consider all k dependencies of W/3k +Ak( x ) to ensure phase transformation invariance 
of our results. In addition, it follows directly from Eqn. flETBI ) that 

— = ~ l $W M /3k( x ) - i r «/3(k) w Q k(x) , (B.7) 
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which is consistent with Eqn.(U), since the fast component of Bloch functions w a k( x ) 
and w^ k (x) are orthogonal in the unit cell. 

Using in Eqn. (|B.5|) the derivative of the expression fl5|), done with the help of 
Eqn. (|B.7|) , gives the exact result 

Im [ i 7^(k,k') p^k, k" 



1 7/3/3 (k, k' 

P^;/3/3( k > k ') = E { r pM 7«/3(k, k') - 7j9a (k, k') r^(k' 



(B.8) 



a^3 



Eqn.( |B.6| ) may now be used to obtain the kp-perturbation expression for r / g a (k) 



r/3a(k) 



X 



3 i(^(ko)-€ QQ (k ))-Ak 
r/3a(k ) + ( Ak • Wko)) r cra (ko) 

77^/3 o-^a 



(B.9) 



Similarly can be obtained the expression for 7^ (k, k') 



7/3«(k,k' 



i(^(ko)-Ak-e aa (ko)-Ak ) 



idpa - i Ak • r <m(k ) <W + i Ak' ■ r /3<x( k o) 5™ (B.10) 



where Ak = k — ko, Ak = k — ko. Upon inserting these results in the expression 
JjD and making use of the vector identity ax (b x c) = b (a • c) — c (a • b), we 



finally arrive at the linearized expression for the scattering shift vector 

R s;/3/3 (k, k') « 
where the vector fig(k ) is defined as 



(k - k ) x ty^ko) , 



^(ko) = V fc x ^(ko) = i J2 r /3a( k o) x r a/ 3(k ) 



(BAD 



(BA2) 



The last expression for x ^(ko) is also derived in Ref.pOfl (Eqn.13). A standard 
analysis using the symmetry properties [H], [T7] of the r a/3 (k ) shows that ng(k ) 
represents an axial pseudo- vector [ID| and, as a consequence, can be nonzero only if 
the material does not posses a center of inversion. It can be evaluated numerically 
from ab-initio values of the respective matrix elements. 
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